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Abstract. We propose a new method to construct an isotropic cel- 
lular automaton corresponding to a reaction-diffusion equation. The 
method consists of replacing the diffusion term and the reaction term 
of the reaction-diffusion equation with a random walk of microscopic 
particles and a discrete vector field which defines the time evolution of 
the particles. The cellular automaton thus obtained can retain isotropy 
and therefore reproduces the patterns found in the numerical solutions 
of the reaction-diffusion equation. As a specific example, we apply the 
method to the Belousov-Zhabotinsky reaction in excitable media. 



1. INTRODUCTION 

Reaction-diffusion equations are used extensiveiy for modeling pattern 
formation observed in natural and social phenomena. The equations are 
deduced from the simple idea that concentration changes of a material in a 
system are caused by reactions between the materials contained in the sys- 
tem and by diffusion of each material. Numerical solutions of the equations 
show a variety of patterns, and they are used widely in various fields such 
as chemistry, biology, medical science etc. pQ. 

An alternative approach to modeling pattern formation is to use a cellular 
automaton (hereafter abbreviated as CA) [2]. CAs are mathematical models 
with discrete time, space and state variables, and are defined by simple time 
evolution rules. They are able to reproduce complex patterns, and therefore 
have a lot of applications. For example, the lattice gas CAs represent very 
well the various features of fluid dynamics and reaction-diffusion phenomena 
[3] . In general, CAs have the advantage of being able to simulate phenomena 
at lower computational cost than when using partial differential equations. 

A lot of models for cooperative phenomena of excitable media, derived by 
using partial differential equations or cellular automata have been proposed. 
The so-called " oregonator" model [SJ EJ [7] expressed in terms of differential 
equations can explain the Belousov-Zhabotinsky(BZ) reaction [4] known as 
a typical example of a reaction-diffusion system. On the other hand, some 
approaches by means of cellular automata for excitable media have been 
proposed in [S] |U] QJJ]. However, these early CA models featured update 
rules based on nearest neighbor connections, and hence, faced several seri- 
ous shortcomings [IT] . The most serious of these is the lack of curvature 
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and dispersion effects and unwanted anisotropy of the front motion |12l [T3] . 
To overcome these problems, several automata models have been proposed 
[13 US EES EGU El EH!- Gerhardt, Shuster and Tyson introduced bigger 
neighborhoods to model the curvature effects and make the threshold a lin- 
ear function to take dispersion into account |12 | ll4 t [To]. Weimar, Tyson and 
Watson improved the models by introducing a mask, i.e. a weighted summa- 
tion of automaton values over large neighborhoods \17\ [T8] , and Henze and 
Tyson extended it to three spatial dimensions |19j . These models recover 
curvature and dispersion effects well, but the anisotropy of wave propagation 
is not completely eliminated. 

In cellular automata, to recover the isotropy of the time-evolution patterns 
is in fact a difficult problem. If the adopted lattice used in the modeling 
has periodicity, such as a square lattice or a hexagonal lattice, then, due to 
this periodicity the time-evolution patterns obtained from the model become 
anisotropic. In |16j Markus and Hess have proposed an isotropic model for 
excitable media. In their model, they use a square lattice, but instead of 
placing each grid point at the center of a unit cell, they assign each grid 
point to a random location within its unit cell. The isotropy is recovered by 
taking a large number of neighboring cells within a circular area with radius 
R. As other CA models for excitable media which seem to recover anisotropy 
to some extent, one can cite the "Moving Average CA" method by Weimar 
[21], the lattice gas method CAs by Raymond et al. [22] or Chen [23] et al., 
the isotropic CA model for the growth process of a bacterial colony by using 
a Voronoi lattice |24] proposed in |25j . In the previous paper |20j . Tanaka 
and the authors also considered an isotropic CA model for the BZ reaction. 
This simple square lattice model adopts a Moore neighborhood with radius 
1 as for the neighborhood cells and recovers the isotropy by introducing 
spatial randomness in relation to a threshold controlling the excitation of 
reaction. 

In the present paper, we propose a new method to construct a CA model 
corresponding to a reaction-diffusion system. The method is based on the 
random walk and on the phase diagram of reaction equations. We require 
the method to satisfy the following two conditions: (i)The time evolution 
pattern of the CA preserves the isotropy of the reaction-diffusion equation. 
(ii)The CA model explicitly contains the control parameters corresponding 
to the reaction terms and diffusion terms in the reaction-diffusion equation. 

In the next section, we explain the general method to construct the 
CA corresponding to a given reaction-diffusion equation. In section 3, the 
method introduced in section 2 is applied to the BZ reaction, and the ob- 
tained evolution patterns are examined. A summary or the results is given 
in section 4. 
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2. METHODOLOGY 

In this section, a method to construct the CA corresponding to a given 
reaction-diffusion equation is introduced. The idea is very simple: to replace 
the diffusion term by a random walk process and the reaction terms by time 
evolution along a discrete vector field obtained from the phase diagram. 
First, we explain the reaction-diffusion equation in brief, then we introduce 
our CA model. 

2.1. Reaction-Diffusion Equation. Suppose that there are N different 
reactive materials U±, U2, • • • , f/jv in some spatial region. Let the densities 
of these materials at position r and time t be Ui(r, t), 112(1", t), • ■ • , ujv(r, t) 
respectively, and put u :=(«i, «2> ■ • • > un) T € K • The reaction-diffusion 
equation is given as 

(i) | = /W + wV 

where D = diag(di, d,2, ■ • • , djv) i.e. D is a N xN diagonal matrix which has 
the diffusion coefficient dj for each material as a diagonal element. The vec- 
tor f(u) = (f\(u), f2(u), • • • , /at(m)) t expresses the interactions between 
the materials and V := is the nabla symbol, in particular V 2 = + ^2 
in two spatial dimension: r = (x,y). 



2.2. The CA Model. We now introduce our cellular automaton model 
for the reaction-diffusion equation ([1]). For simplicity, we assume that the 
equation is defined in two spatial dimensions and that the CA model is 
defined on a two dimensional square lattice. Generalization to higher di- 
mensional systems and different types of lattices will be straightforward. 
In the reaction-diffusion equation ([T|), the variable u represents the density 
of reactive material. In our CA model, we replace the density of material 
u(r,t) by the number of microscopic particles it^ mn € 1> + N at the corre- 
sponding lattice point (m,n) E 1? at time step t 6 Z/2. Then the time 
evolution of our CA model is determined by 



u t+l/2 i , R{ii t N 

u t+1 = u t+l/2 + F(u t+l/2 ) 

where u t mn denotes a set of concentration variables around (m, n) at time 
step t, i?('uj rm ) and i^it^) denote discrete vector fields which are obtained 
by replacing the diffusion and reaction terms of the reaction-diffusion equa- 
tion respectively with the following processes: 

(i) The diffusion term L>V 2 w corresponds to a random walk of particles. 
Let us consider a random walk of particles defined in the Neumann neigh- 
borhood u l mn = {u l m n , u l m _ hn , w' m+ljn , }• One can equally 
adopt the Moor neighborhood or other neighborhoods. We denote by U^ n 
the stochastic variable which defines the number of particles moving from 
the (m, n) site to the (m + 1, n) site due to the random walk of the particles, 
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by U^t ln that from from (m — l,n) to (m, n), and so on. Then i?('U^ m ) 
equals the difference between the number of outgoing particles and that of 
the incoming particles due to the random walk at position (m, n) and time 
t: 

R( u mn) = ^m-l,n + ^m+l,n + ^m,n—l + ^m,n+l 
V u / V*-' mn 1 ^ mn 1 ^ mn 1 ^ mn)' 

If we define pj as the transition probability of particles to one neighbor- 
ing cell in the random walk, 1 — 4pj is the probability of particles to stay 
on site and the expectation of U^ n , U^-i n e ^ c - are given respectively 
by = Pu l m _ X n , (J7^i ljn ) = Pu l m _ l n with a diagonal matrix P = 

diag(pi,P2i " " ' )Pn)- We can see that P corresponds to the diffusion coeffi- 
cient D of the reaction-diffusion equation. 

(ii) The reaction term f(u) is replaced with an appropriate discrete func- 
tion = (FiO^J, F 2 (u t mn ), • • • , F N (u t mn )) T € Z N . In equation ([I]), 
the reaction term f(u) is the vector field that defines the velocity vector 
Hence F(x) should be so chosen such that the time evolution of x is consis- 
tent with the typical orbits in the phase diagram of the ordinary differential 
equation 

du 

!Tt= f{u) - 

Since F(vr mn ) sometimes returns a negative number, it is practically con- 
venient to use the discrete vector field, G(tt^ n ) := u mn + F(v!" mv ) > 0. 
Then our CA model is rewritten as 

u t+i/2 =u t ■ R(il t \ 

"ran V "mn / • 



We can obtain the time evolution pattern by successive substitution of 
Riulnn) and G(u^ n ) for appropriate initial conditions. 

It is well known that the continuous limit of the random walk is equiva- 
lent to a diffusion equation, and that in the large scale limit, isotropy of the 
distribution function of the particles is guaranteed by a random walk. Since 
the discrete vector field F(u) is chosen such that it is essentially equivalent 
to the vector field f(u) in the continuous limit, we expect that the patterns 
obtained from our CA model become almost isotropic in certain large sys- 
tems. It should be noted that our model naturally contains the parameters 
{Pi} which correspond to the diffusion coefficients and all other control pa- 
rameters for reactions, which are necessarily contained in the discrete vector 
field. 



3. CA MODEL FOR THE BeLOUSOV-ZhABOTINSKY REACTION 



In this section, we apply the method introduced in the previous section 
to the Belousov-Zhabotinsky (BZ) reaction as a specific example. First, 
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we briefly explain the BZ reaction and oregonator known as a mathematical 
model for this reaction. Next, we introduce our CA model of the BZ reaction. 

3.1. BZ Reaction. The BZ reaction is known as an oscillating oxidation- 
reduction reaction which occurs by mixing some chemical compounds (such 
as Ce 4+ , BrO 3 ", CH 2 (COOH) 2 , H 2 S0 4 ). If the BZ reaction is spread spa- 
tially, then it forms trigger waves, spiral waves or target patterns The 
BZ reaction is often modeled by using partial differential equations. Among 
them, the oregonator, which is a system of simultaneous ordinary differen- 
tial equations with two variables, is widely considered to be the simplest 
possible model [6]. The time evolution of the spatial patterns in the BZ 
reaction is described by the following equation which adds diffusion terms 
to the oregonator: 

du 
di = 



bv(u-a) 

u 1 — u) 



u + a 



+ d u V 2 u, 



(5) 

OV . _n 

— = u - v + d v V v, 
ot 

where b and e (or 1/e) are, respectively, a threshold which gives the excita- 
tion and a parameter which defines the excitability of reaction respectively. 
Depending on the parameters a, b and e, Eq. © shows two typical states: 
an excitable state with one stable equilibrium point and an oscillatory state 
with one unstable equilibrium point. Here we consider only the excitable 
state (or excitable media). 

Let f(u,v) := ^ [u(l — u) — ^^_ ] and g(u, v) := u — v. Figure [JJ shows 
the phase diagram for the excitable state of the BZ reaction. The null 
clines that are obtained from f(u,v) = g(u,v) = are shown by solid lines 
and a typical orbit for the excitable state is shown by a dashed line. The 
intersecting point of / = and g = is a stable point if there is no diffusion, 
however it becomes unstable when the strength of the perturbation (mainly 
due to the diffusion effects) exceeds a certain threshold S. Then the state 
of the medium becomes unstable and changes along the dashed line shown 
in the phase diagram, until it returns to the equilibrium point once again. 
The repetition of this process induces spacial patterns such as spiral waves. 

3.2. The CA Model. Our CA model for Eq. © is described in the form 
of Eq. (jl]), introduced in section 2.2. According to Eq. ([5]), we put 

u mn := { u mni v mn) ) G(ll mn ) := (G u (u mn ) , G v {u mn )) and R(u mn ) '■ = 

(R u (u mn ), R v (v mn )) . 

First, we define the discrete vector field G(u^ n ) by imitating the solution 
orbit of the phase diagram shown in Figure [TJ The discrete vector field we 
configured is shown in Table [T] and illustrated in Figure [21 We can see that 
Figure [2] is a simplification of the phase diagram of Figure [TJ Let u t mn G Z + , 
v mn ^ {0' 1} i n tlfi s example of an excitable BZ reaction. Parameters a, 
(3 and 7 control the rate of reaction, and A is the threshold which deter- 
mines whether the excitation occurs or not. Parameter iV £ denotes the 
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u 



Figure 1. The phase diagram for the excitable BZ reaction. 

expected maximum value of variable u t mn . All of these five parameters are 
positive integers. 

The value of the discrete vector field G(w^j n , u^ n ) is determined in func- 
tion of the values of u t mn and v t mn , as shown in Table [TJ For < u t mn < 
A, v t mn = 0, the state of the medium returns to the equilibrium point with 
velocity a, because the state cannot exceed the threshold A due to inade- 
quate diffusion effects. For A < -u^ n < N — 1 — /3, v t mn = 0, the variable u t mn 
increases at the rate j3 until u t mn > N — 1— f3. In the range 7 < !!*„, v t mn = 1, 
the variable u t mn decreases at the rate of 7. In the two remaining ranges, 
the variable v l mn changes between and 1. Here we assumed v t mn € {0, 1}, 
because two states are sufficient to separate the phases, according to the 
medium, corresponding to the variable u^ n . 



Table 1. Example of discrete vector field Gr( , uf rm ) for an 
excitable medium. 



("mni "mil) ^mn 


(G u (u mn ),G v (u mn )) — G(u mn ) 


(0<n s mn < A,v t mn = 0) 


(max[u^j n — a, 0], 0) 


(A iV- 1 = 0) 


(4 + P, 0) 


(N-l-f3<u t mn ,v t mn = 0) 


(iV-1,1) 


(7 *C u m ni ^mn -Q 


( u mn ~ 1i 1) 


(0 < u l mn < 7,< n = 1) 


(0,0) 



In Eq. (JSJ) , the excitation of excitable media is determined by the diffusion 
of u, i.e., d u ^> d v . Hence, in this example, we only consider the random 
walk of the variable u l mn and we simulate the time evolution by putting 
P = diag(p u ,0). 
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FIGURE 2. Outline of the discrete vector field G{u mn ) in TableHJ 
Consequently, our CA model for the excitable BZ reaction is given as: 



= M t 

(6) 



u mn' u mn ^ui^rn 

v t+1/2 = v l 

u mn ^mm 



,.t+l _ Q ( u t+l/2\ 

v t+l _ q ( U t+V2) 
u mn ~ ^vyu-mn )• 

3.3. Numerical Results. In this section, we show several time evolution 
patterns obtained from Eq. © and discuss the results. 

3.3.1. Patterns. The first example is a single trigger wave (Figure [3]). It 
is produced by the initial condition u^ nn = [h ■ exp(— ((to — L/2) 2 + (n — 
L/2) 2 )/w 2 )] and v^ nn = on a 2-dimensional square area with L x L cells. 
Here [ ] is Gauss' symbol, i.e.: [x] is the largest integer that is less than or 
equal to x and h and w are positive real numbers. We observed that from 
a pulse triggered at the center, a ring-shaped wave spreads outwards in an 
almost isotropic fashion. 

In order to generate a spiral wave, appropriate initial conditions are nec- 
essary. Firstly, we generate a single trigger wave like the one shown in Figure 
[31 Then we cut off one part of the ring pattern as shown in Figure H] and 
use the remainder as the initial state. The spiral wave obtained from our 
model is shown in Figure [SJ 

The third example is a target pattern (Figure [6]). The initial condition is 
the same as in the case of the trigger wave, but with different parameters. 
From the central pulse, ring-shaped waves are produced repeatedly. 

3.3.2. Anisotropy. We evaluate the anisotropy of the trigger wave in Figure 
[3] by measuring the residual error when compared with the average radius of 
the ring, and plot it in Figure [7] as a function of the propagation direction. 
We find that the wave fronts of the trigger wave propagate in each direction 
with an anisotropy in the range of ±2.5 percent. In Figure |H] we plot the 
variation from the circle as a function of the radius of the trigger wave. It 
shows that the trigger wave indeed grows closer to a complete circle as the 
radius (and therefore time) increases. 
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Figure 3. Single trigger wave. 500 x 500 cells, t = 5850, 
N = 100, p u = 0.2, A = 21, a = 1, (3 = 1, 7 = 1. 




Figure 4. A cut trigger wave. 200 x 200 cells, t = 250, 
N = 30, p u = 0.2, A = 6, a = 1, f3 = 2, 7 = 1. 

Next, we evaluate the parameter dependency of the anisotropy for the 
patterns observed from our model. In Figures [9] and [TUl the variations of 
the wavefronts of trigger waves propagated to radius 450, are plotted as a 
function of the transition probability p u of particles by changing parameters 
P (Fig. HJ) and A (Fig. [101) . We find that the patterns become more isotropic 
when diffusion becomes stronger, and that they tend to be isotropic for 
smaller j3, but depend less on A. 

3.3.3. Parameters. Our CA model has six parameters N, p u , A, a, f3 and 
7. The parameter A corresponds to b in Eq. ([5]) and defines the threshold 
of excitation; j3 corresponds to 1 /e which defines the excitability of reaction. 
Figure [TT] is the phase diagram of our model obtained by changing the values 
of A and f3. The remaining parameters are chosen as iV = 50, p u = 0.1, a = 
7 = 1. We confirm that spiral waves can be generated in a wide region of 
the parameters. In an even larger parameter range, we find there exist three 
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Figure 5. Spiral wave. 200 x 200 cells, t = 963, N = 30, 
p u = 0.2, A = 6, o = l, /3 = 2, 7 = 1. 




FIGURE 6. Target pattern. 300x300 cells, i = 626, N = 100, 
p u = 0.04, A = 2, a = 1, p = 10, 7 = 1. 

different regions; one without propagating waves, one allowing for trigger 
waves (or broken trigger waves), and one chaotic pattern region. 

The parameter range for the spiral waves obtained from Eq. ([5]) has been 
examined in detail by Jahnke and Winfree [7j. The general tendency is that, 
when b is sufficiently large, the spiral wave does not appear and for l/e> 1 
the spiral wave destabilizes and chaotic behaviour is observed. As shown in 
Figure [Til the phase diagram of our model has the same feature as that for 
Eq. (J3J), and we may conclude that the parameters A and /3 play the same 
role as the control parameters of the reaction terms in the reaction-diffusion 
equation. 

4. SUMMARY 

We have proposed a method to construct a CA model corresponding to 
a reaction-diffusion equation, in which the diffusion effect is replaced by 
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Figure 7. Plot of the anisotropy of the model, measured as 
the residual error of the trigger pattern in Figured 
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Figure 8. The relationship between the radius of the ring 
and the variation of the radius. 

a random walk with transition probability matrix P and the reaction by 
discrete vector fields. The model can include control parameters for both 
diffusion and reaction, as is the case in the reaction-diffusion equation. As 
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Figure 9. The parametric dependency of the anisotropy of 
the model for the parameter (3. N = 100, A = 15, a = 1, 
7 = 1. 

an example, we have shown that our model can successfully reproduce the 
patterns of BZ reaction. Applications to other reaction-diffusion systems 
such as FitzHugh-Nagumo equation |26| are interesting future problems. In 
the present method, however, there are still various candidates for the time 
evolution rules, depending on the choice of the discrete vector fields that are 
supposed to have similar features to those of the continuous vector fields, 
given by the reaction-diffusion equation. Of course we can also adopt dis- 
crete vector fields by investigating the reaction process from a microscopic 
point of view. The determination a suitable time evolution rule will de- 
pend on the individual phenomenon and we will investigate this problem 
extensively in the future. 
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Figure 10. The parameter dependence of the anisotropy of 
the model on the parameter A. N = 100, a = 1, /3 = 1, 
7 = 1. 
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